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Abstract 


A two-time scale asymptotic method has been introduced to analyze the mul¬ 
timodal mean-field Kuramoto-Sakaguchi model of oscillator synchronization 
in the high-frequency limit. The method allows to uncouple the probabil¬ 
ity density in different components corresponding to the different peaks of 
the oscillator frequency distribution. Each component evolves toward a sta¬ 
tionary state in a comoving frame and the overall order parameter can be 
reconstructed by combining them. Synchronized phases are a combination 
of traveling waves and incoherent solutions depending on parameter values. 
Our results agree very well with direct numerical simulations of the nonlinear 
Fokker-Planck equation for the probability density. Numerical results have 
been obtained by finite differences and a spectral method in the particular case 
of bimodal (symmetric and asymmetric) frequency distribution with or with¬ 
out external field. We also recover in a very easy and intuitive way the only 
other known analytical results: those corresponding to reflection-symmetric 
bimodal frequency distributions near bifurcation points. 
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I. INTRODUCTION 


In recent years mathematical modeling and analysis of synchronization phenomena re¬ 
ceived increased attention because of its occurrence in quite different helds, such as solid 
state physics M. biological systems [§-0, chemical reactions pj, etc. These phenomena 
can be modeled in terms of populations of interacting, nonlinearly coupled oscillators as 
hrst proposed by Winfree [Q. While the dynamic behavior of a small number of oscilla¬ 
tors can be quite interesting |^, here we are concerned with synchronization as a collective 
phenomenon for large populations of interacting oscillators [^. Then we can describe pop¬ 
ulations of oscillators interacting via simple couplings (e.g., all-to-all, mean-held couplings) 

isiini 


by means of kinetic equations for one-oscillator densities 


Many recent studies of 


synchronization phenomena combine numerical simulations with linear stability and bifurca¬ 
tion analyses of particular (stable) incoherent and synchronized states IM- These works 
have described the onset of synchronized phases and, near degenerate bifurcation points, 
synchronized phases from their beginning to their end in the corresponding bifurcation dia¬ 
gram |]I^. In this paper we introduce a high-frequency singular perturbation method which 
describes (in a conveniently analytical manner) synchronized phases far from bifurcation 
points. The method nicely agrees with the results of numerical simulations. 

These ideas may be made concrete in a simple model put forth by Kuramoto and Sak- 
aguchi (see also P). It consists of a population of coupled phase oscillators, Oj(t), 

having natural frequencies ujj distributed with a given probability density g{uj), and subject 
to the action of an external held hj which is distributed with a probability density f{h), 


N 


9j = ojj -t- ij{t) — hj sinOj + ^ Kji sm{9i — 9j), 




1=1 


Here are independent white noise processes with expected values 

im) = 0, mMt')) = 2D5it - t') 5,1. 


( 1 ) 


( 2 ) 


In the absence of external held and white noise, each oscillator tries to run independently 
at its own frequency while the coupling tends to synchronize it to all the others. When 
the coupling is sufficiently weak, the oscillators run incoherently whereas beyond a certain 
threshold collective synchronization appears spontaneously. So far, several particular pre¬ 
scriptions for the matrix Kji have been considered. For instance, Kji = K > 0 only when 
\j — l\ = 1, and Kji = 0 otherwise (next-neighbor coupling) [|^; Kji = K/N > 0 (mean-held 
coupling) |1l6|,^; hierarchical coupling |^; random long-range coupling |p0|-p^ or even state 
dependent interactions |^. In the mean-held case, the model (|^)-(Q) can be written in a 
convenient form by dehning the (complex-valued) order-parameter 


re 


1 ^ 

piSi 

N ^ 

i=i 


( 3 ) 


Here |r(t)| measures the phase coherence of the oscillators, and ipit) measures the average 
phase. Then eq. (1 reads 


9j = ojj — hj sin 9j + Kr sin('0 — 9j) -\- ^j{t), j = 1,2,..., N. 


( 4 ) 
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In the limit of infinitely many oscillators, iV —> cx), a nonlinear Fokker-Planck equation 
(NLFPE) was derived |T0|JTl|] for the one-oscillator probability density, p{9,t,uj, h), 


dp _ _ d , . 

dt ^ de^ 

the drift-term being given by 

n(6*, t, a;) = ui — hsinO + Krsiniplj — 9), 
and the order-parameter amplitnde by 


re^P = 


r2n r+oo f-\-oo 


e^^p{9, t, u!, h) g{uj) f{h) d9duj dh. 


(5) 

( 6 ) 

(7) 


' —OO 8 —OO 


The probability density is reqnired to be 27r-periodic as a fnnction of 9 and normalized 
according to 


p 2 TT 

/ p{9^ t, uj, h) d9 = 1. 

Jo 


( 8 ) 


Mean-held models such as those described above were studied, e.g., by Strogatz and 
Mirollo 0 in the absence of external held and for a unimodal [g{u}) is non-increasing for 
a; > 0] frequency distribution, g{uj), having rehection symmetry, g{—uj) = g{uj). In |]^, the 


authors showed that for K smaller than a certain value Kc, the incoherent equiprobability 
distribution, po = l/(27r), is linearly stable, and linearly unstable for K > Kc. As D ^ 0-I-, 
the incoherence solution is still unstable for K > Kc [= 2/7rp(0) at D = 0], but it is neu¬ 
trally stable for K < Kp. the whole spectrum of the equation linearized about po collapses 
to the imaginary axis. In [^], the non/mear stability issue was addressed, and the case of a 
rehection-symmetric bimodal frequency distribution was considered [g{oj) is even and it has 
maxima at a; = icuo]. In this case, new bifurcations appear, and bifurcating synchronized 
states have been asymptotically constructed in the neighborhood of the bifurcation values of 
the coupling strength. The nonlinear stability properties of such solutions were also studied 
for the explicit discrete example g{uj) = ^[5{u! — uq) -|- (5(a; -|- cuo)], cf. 


T2| . A complete 


bifurcation study taking into account the reflection symmetry of g{u}) was carried out by 
Crawford, |^. Similar results were obtained by Okuda and Kuramoto in the related case 


of mutual entrainment between populations of coupled oscillators with different frequencies 
0. Furthermore, a two-parameter bifurcation analysis near the tricritical point (at which 


bifurcating stationary and oscillatory solution branches coalesce) allows ns to visualize a 
global bifurcation diagram in which oscillatory solution branches may be calculated analyt¬ 
ically from their onset to their end [^]. The effect of an external held on Kuramoto models 
has been analyzed in Refs. p^,p5 


In this paper we shall illustrate our high-frequency perturbation method by applying it 
to the generalized mean-held Kuramoto model (|^)-(|^). We shall assume that the frequency 
distribution is multimodal in the high-frequency limit: p(a;) has m maxima located at 
/ = 1,..., m, where ujq —> cx). Then p(a;) duj tends to the limit distribution 


F(f 2 ) dn = J2ai 5{n - ni) dn, 


( 9 ) 


1=1 


00 


with = 1, and = 

1=1 ^0 
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independently of the shape of g{uj) as ujq —>• oo. Then g{uj) duj and r(f2) dVL may be used 
interchangeably when calculating any moment of the probability density [including of course 
the all-important order parameter Thus any frequency distribution is equivalent to a 
discrete multimodal distribution in the high-frequency limit. The discrete symmetric bi- 
modal distribution considered in |[T^JT5[1 corresponds to m = 2, ff; = (—1)^ ai = = 1,2. 

We shall show that the oscillator probability density splits into m components, each con¬ 
tributing a wave rotating with frequency VtioJo to the order parameter. The envelope of each 
component evolves to a stationary state as the time elapses. Thus our method yields ana¬ 
lytical expressions for the probability density and the order parameter during the transients 
toward synchronized (or incoherent) phases, which agree with direct numerical simulations 
of the NLFPE. Since it is not a small-amplitude expansion, our method is valid well inside 
the regions of stable synchronized phases in the phase diagram, far from bifurcation points. 
Of course we have derived the method in the limit wq —> cxo, but comparison with numerical 
simulations shows that wq = 7 is already close to inhnity for all practical purposes. 

Our numerical calculations have been carried out by means of hnite difference schemes 
and by using a spectral method which generates a hierarchy of ordinary differential equations 
for moments of the probability density which include the order parameter. This method is 
equivalent to an expansion of the probability density in a Fourier series and it could in prin¬ 
ciple be used to reconstruct it. The moment hierarchy was derived directly from Eqs. ®-(@) 
by Perez Vicente and Ritort [^. They assumed that arithmetic means converged to en¬ 
semble averages in the limit N ^ oo [keeping t = 0(1)], which was justihed in |^. From 
the moment hierarchy, Perez Vicente and Ritort also derived a nonlinear kinetic equa¬ 
tion for a moment-generating function T(6*, y,t) = f p{9, t; a;) g{uj)duj, which is related to 
functional-equation formulations of Equilibrium Statistical Mechanics [^ and fluid turbu¬ 
lence [P8| ,P5[[. 

The rest of the paper consists of a description of our method of multiple scales in Section 
II, comparison with numerical results in Section III, and our conclusions in Section IV. 


II. METHOD OF MULTIPLE SCALES 


The high frequency limit of the NLFPE can be analyzed by means of a method of multiple 
scales. Let us change variables to a comoving frame and therefore rewrite the equation as 

dp d'^p dU p 

m 


dl3‘^ djd ’ 

U = Kr sin('0 — fd — out) — h sm{(d + ujt) 

= Kr sin — (3 — — h sin ^(3 + , 

n 

(3 = 0 — uit = 9 -1, 

£ 

1 

£ = — -Cl. 

UIq 


The order parameter is now 


re 


_ 


i=i 


p(/3, t; Oj, h; e) f{h) d(3 dh, 


( 10 ) 


( 11 ) 

( 12 ) 

(13) 

(14) 
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where we have used and have implicitly assumed that h=0(l). The discrete character 
of the frequency distribution in the high-frequency limit makes it possible to simplify (0. 
In fact, Eq. ( p!l[) shows that p may be split in different components pj = p{P,t-,flj,h;e). 
Therefore we can write (pTl) as a coupled system of equations for the density components pj: 


dpj ^ ^Pi _ djUjPj) 

dt d(3^ d(3 ’ 

Uj = lm\Kj2 ^ 


eh/3'-/3) f. h- 0clh _ h eh/3+n,|) 


1 = 1 


p2'k 

/ Pj(/9, t] h] e) dj3 = 1. 

0 


(15) 

(16) 
(17) 


Eqs. (|T^ and (pT]) contain terms with rapidly-varying coefficients. It is then to be 
expected that an appropriate asymptotic method will be able to average them out thereby 
capturing the slow evolution of pj (or perhaps its envelope). This may be achieved by 
introducing fast and slow time scales as follows: 


t 

£ 


t = t. 


(18) 


We look for a distribution function which is a 27r-periodic function of (3 according to the 
Ansatz: 


p(/3, t; 12, h; e) = p^^\(3, r, t; 12, h) + ep^^\l3, r, 2; 12, h) + 0{e^), 
Inserting (P!U|) into (HD-®, we obtain the following hierarchy of equations 


where 


Eq. 




5r = °- 


(19) 


( 20 ) 


dp. 


( 1 ) 


dr 


^ \ pf Im I a: ^ aj - h 


d(3 

dpT 


= + - Kc, A {pf Im (e-«>zf)} , 


dt 


df3^ 


d(3 


Zf\t) = / p'f\p,t,h) f{h) dpdh. 


do) 


implies that p^p is independent of r. Then the terms in the right side of (^T|) 
which do not have r-dependent coefficients give rise to secular terms (unbounded on the 
r-time scale). The condition that no secular terms should appear is 


( 21 ) 


( 22 ) 


dp 


(0) 


^ + Ka, A {pf I„, »>)} = 0. 


dt df3^ ^ df3 
This equation should be solved for pf'^ together with (1^) , the normalization condition 


(23) 
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( 24 ) 


f pf\P,t-,h)dl3=l, 


and an appropriate initial condition. We see that, except for the h-integration in this 
problem is eqnivalent to solving a NLFPE with frequency distribution g{uj) = (identical 
oscillators) and coupling constant Kj = Kaj. If the initial condition is independent of 
the external held h, we know that the solution of the previous NLFPE evolves towards a 
stationary state as time elapses [^. If the initial condition depends on h, all we can say is 
that / h) f{h) dh tends to a stationary state independent of h as t —>• cxd. In both 
cases all possible stationary states are solutions of Eqs. (^^-(^) below 

pKajRjD~^ cos{^j—0) j2Tr ^—KajRjD~^ cos(^j—/3—/3') 


pf(0) = 


2-k 


icos('i'j /3) e 


cos('I'j—/3—/3') 


(26) 


The order parameter is calculated by inserting (^5[) into (^): 




p2-K 


^ (rj) dr] = lim Z. 

t^OC 


(0) 


it). 


(26) 


For Kaj < 2D, the only stationary solution is ^ (incoherence), which is stable. 

At Kaj = 2D, a stable branch of synchronized solutions bifurcates supercritically from 
incoherence. They exist for all Kaj > 2D. 

The overall order parameter (|T^ is given by 


rDi’ = ^ aj Rj + 0{e ). 

i=i 


(27) 


To hnd 'ip, we multiply both sides of (^Tj) by e and then take imaginary and real parts. 
After a little algebra, we obtain 


sin(fi,r + T,) 

ET=i Rj cos(^r + T,) ’ 

(28) 

r = '^aj Rj cos{fljT + Tj — jp). 

(29) 




Notice that r in (^^ may be negative, positive or zero. Then the amplitude of the overall 
order parameter is |r(t)|. 

Let us now consider, for the sake of dehniteness, the special case of an asymmetric 
bimodal frequency distribution, with zero external held. 


F(o;) = a 5(f2 - 1) + (1 - a) 5(fi + 1), 0 < a < 1, /(h) = 5{h), (30) 

and analyze the possible synchronized states. Eqs. (PSD and ( PD|) become 

, _ sin(T++ r) + (1 - a) sin(T„ - r) 

a R+ cos(T+ + r) + (1 - a) R. cos(T_ - r) ’ ^ ^ 

r = aR^ cos(T+ + r — />) + (1 — a)i?_ cos(T_ — t — ip). (32) 


Let us now assume that a < 1/2 to be specihc. Then we have the following possibilities 
depending on the value of the coupling constant: 


6 








1 . IfO < K < 2D/(I — a), the incoherent solution p = l/(27r) is stable and it is the only 
possible stationary solution. 

2. If 2D/{I — a) < K < 2D/a, a globally stable partially synchronized solution issues 

forth from incoherence ai K = 2D/{I — a). It has = 0, 'ip = — r, and 

r = (1 — a)R-. Its component p+ = l/(27r) is incoherent while its component p_ is 
synchronized according to Eq. The overall effect is having a traveling wave solu¬ 
tion (rotating clockwisely), once the angular variable (3 is changed back to 9 according 

to(|T2D. 


3. If > 2D/a, the component becomes partially synchronized too. The probability 
density then has traveling wave components rotating clockwisely and anticlockwisely. 
Their order parameters have different strengths and R- > R+ if a < 1/2. 


When a = 1/2, both traveling wave components appear at the same value of the coupling 
constant, K = AD, and have equal strength: R^ = R_ = R, T_|_ = = T. i? is the 

amplitude of the order parameter corresponding to a unimodal frequency distribution and a 
coupling constant K+ = = aK = K/2. Then (|3l[) and imply that V’ = ^ + <? 7 r (g is 

an integer number) and r = R cos{ujot + qn) , respectively. Thus we have obtained an overall 
standing wave which is stable. Of course other possible solutions are traveling waves with 
> 0 , = 0 and = 0 , > 0 , which should be unstable because incoherence is an 

unstable solution of (pS] ) for the corresponding stationary component pj if K/2 > 2D. These 
results coincide perfectly with those obtained by means of bifurcation theory in and 
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for a symmetric bimodal frequency distribution (a = 1/2). To see this, we recall that the 
stable (up to a constant shift in the origin of time which depends on initial conditions) 
standing wave probability density may be approximated near a bifurcation point = AD 
by the following expressions ||T^: 


p{6,t,u}) “ 


K = AD+ e^ K2, 


ai = A 




p{Ut—9) 


+ 


( H-|-i(f2-|-cu) D + i{D — (j) 
A = 


Re a 


iue^K2t 


+ cc, VL = \Jujq — D‘^, 
Im (7 - 1 - P) 


Re (7 + p) 


V = Im a 


Re (7 + p) 


Re a , 


(33) 

(34) 

(35) 


where cc means complex conjugate of the preceding term. As cuq —»• +cxo, the parameters a, 
P, 7 of (|35D become |T5| 




(36) 


Inserting (^3]) to (p6|) in Eq. O for the order parameter, we obtain that p is constant and 


K - AD 

r{t) = W —^ 2 )— cosuat + 0((K - /C)) ■ 


(37) 


Now we can compare Eq. ( 0 ) with the result of our two-time scale method, r = R cosujQt. 
R is the amplitude of the order parameter corresponding to a unimodal frequency distri¬ 
bution and a coupling constant = K_ = aK = A'/2. Near the bifurcation point. 
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Eq. (2.12) of Ref. with Kc = 2D (corresponding to ojq = 0) and K± = 7^/2 yields 
R ~ y {K/2 — 2D)/D, which implies exactly the result (^71). It is immediate to show that 
both methods also lead to the same expressions for traveling wave solutions. 


III. NUMERICAL RESULTS 
A. Spectral numerical method 

Direct numerical simulations of the nonlinear Fokker-Planck system conhrms our asymp¬ 
totic results. We have studied discrete bimodal frequency distributions only, and used two 
different numerical methods. A standard hnite differences method may be used to numeri¬ 
cally integrate (j^) - (H) without stability problems up to frequencies ujo = 15 (we set D = 1 in 
all our computations). For larger frequency values, time steps below 0.008 were needed and 
the computing cost makes this method unpractical. As indicated in the previous Section, 
the drift term dominates diffusion at higher frequency and the system acquires a quasi- 
hyperbolic character. To simulate the NLFPE at high frequencies, we propose a simple 
spectral method, which we will describe in the simple case of f{h) = 6{h). The idea is to 
hnd a set of ordinary differential equations for moments of the probability density related 
to the order parameter 


™(i) _ 
x^ — 

p27T 

/ p{6, f, iwo) cos[j{i> - 6^)] de, 

Jo 

(38) 

y± = 

p27r 

/ p{9, t, icuo) sin[j('0 - 6^)] d9, 

JO 

(39) 


(1) 1 n ^ (1) 

r = axf -1- (1 — a) x_ . 

(40) 


An inhnite hierarchy of equations for these moments may be obtained by differentiating 
(^) and (1^) with respect to time and then using the NLFPE and integration by parts to 
simplify the result. We obtain 


dt 

dy^± 

dt 


-fx‘'£±ju}oy±^ + 


D ..O') 


Kj 


Ky 

(i+1) 


dip 

U) 

(41) 

2 

/y^ /-yy ' 

2 

rxf 

- i 

dt 

y±R 

Kj 

O'-I) 

Ky 

(i+i) 

+ j 

dip 

U) 

(42) 

2 

ryf 

2 

ryf 

dt 

, 


II 

= (Uo [q 

— 

(1- 

a) 


(43) 


As explained in the Introduction, an equivalent hierarchy may be derived directly from the 
Langevin equations (|l]) and (H) 


The numerical method consists of solving 


for 


j = 1,..., A, with = 0. The number of modes, N, should be chosen so 

large that the numerical results for the order parameter do not depend on it. A practical 
case is presented in Fig. ^ for (Uq = 15, iP = 6 and a = 0.5 for which the method of hnite 
differences is still practical. We see that keeping four modes {N = 4) yields already quite 
good agreement. Let us now describe the results of our numerical simulations. 
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B. Results for bimodal frequency distributions and no external field 


We see in Figures ^ to ^ that our analytical (asymptotic) and numerical results agree 
very well except for a constant phase shift which decreases as Uq increases (compare figures ^ 
and ^ corresponding to ujq = 15 and 200, respectively). Results for an asymmetric bimodal 
frequency distribution without external field are depicted in Figs. ^ to p. As explained 
in the previous Section, we obtain different synchronized phases depending on the value 
of the coupling constant for each component of the probability density. In Figures ^ and 
I, K > 2D/a > 2D/{I — a). Then each component of the probability density evolves 
towards a synchronized phase rotating with its own frequency, ±ct;o, and with a constant 
amplitude of the order parameter given by the stationary expression (1^). The overall order 
parameter is given by Eqs. (PTD-(P2[) and the difference between analytical and numerical 
results is a constant shift in time which diminishes as the frequency uq increases. In Fig. ^ we 
observe the situation for a smaller coupling constant such that only one density component is 
synchronized. We obtain a traveling wave whose order parameter has a constant amplitude 
and a phase linearly decreasing with time. What happens if the frequency distribution has 
reflection symmetry {a = 0.5) is obvious: both density components have equal strength and 
therefore the phase of the order parameter is constant and its amplitude oscillates giving rise 
to a standing wave. This is exactly what bifurcation theory predicts |^,^. We have checked 
the excellent agreement between our present asymptotic theory, the leading-order expression 
for the order parameter obtained by bifurcation theory, and direct numerical simulations. 
The results obtained by these three methods are indistinguishable for K = 4.005 {Kc = 4). 


C. Results for unimodal frequency distributions and deterministic external field 


Our asymptotic method yields analytical results when external fields of magnitude small 
compared to ujq are included. For the sake of simplicity we shall present results correspond¬ 
ing to unimodal frequency distributions, g{(jj) = 5(a; — (Uq), and external field distributions, 
f{h) = 6{h — ho). Then the probability density has a unique component rotating at fre¬ 
quency ujq which evolves toward the stationary distribution (^) (in the comoving frame). 
This prediction is qualitatively supported by the numerical simulations as depicted in Fig. |^. 
The numerical results show that the amplitude of the order parameter oscillates about the 
constant value predicted by our asymptotics. The difference is of order £ and it could be ac¬ 
counted by the first correction to the leading-order result. Fig. ^ shows that the oscillation of 
the order parameter amplitude is enhanced as hg increases. Finally all oscillations disappear 
if the external field becomes of the same order as the frequency, as depicted in Fig. |^. Notice 
that our method supports [in the limit (Ug —>■ oo, h = 0(1)] a conjecture by Arenas and Perez 
Vicente [^ 


the amplitude of the order parameter in the oscillatorily synchronized state 


m 


the presence of an external field) is given by the same expression as in the stationary state if 
the exact time-dependent phase of the order parameter is inserted (instead of the stationary 
phase). Of course it seems that the conjecture holds for a wide variety of parameter values, 
some outside the range of validity of our asymptotic method 


25 
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IV. CONCLUSIONS 


The high-frequency limit of the mean-held Kuramoto-Sakaguchi model of oscillator syn¬ 
chronization has been stndied by new multiscale and nnmerical spectral methods. The main 
result of the multitime scale method is that the probability density splits into independent 
components corresponding to the different peaks in the oscillator frequency distribution. 
Each density component evolves towards a stationary distribntion in a comoving frame 
rotating with the frequency of the corresponding peak in the oscillator frequency distribu¬ 
tion. The overall order parameter may be calculated by putting together the partial order 
parameters of different components. This gives a simple picture of overall oscillatory syn¬ 
chronization by studying synchronization of each density component. Onr method gives the 
same results as bifurcation theory for those parameter valnes where both approximations 
hold. Onr asymptotic method also works far from bifnrcation points and it agrees well with 
results of numerical simulations. We have used a spectral method to numerically integrate 
the nonlinear Fokker-Planck system for frequency values where simple hnite differences break 
down. This method has merit in itself and should be studied in more detail separately. 
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FIGURES 


Fig.l 



FIG. 1. Comparison between the results of numerical simulations by finite differences and our 
spectral method. We have a discrete bimodal frequency distribution of the oscillators, no external 
field and the following parameter values: loq = 15, K = 6 and a = 0.5 (frequency distribution with 
reflection symmetry). Differences between the methods are appreciated only on a rather fine time 
scale as the inset shows. 
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of Figure 1: (a) analytical result from our leading-order asymptotics; (b) numerical simulation; (c) 
comparison between both results; (d) same as (c) but now we have shifted the analytical result so 
that t —> t -|- 0.055. 
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FIG. 3. Same as Fig. 2(c) except that now loq = 200. Notice that the time shift between 
analytical and numerical results is now much smaller. 
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(d) 



FIG. 4. Time evolution of the order parameter r{t) for an asymmetric bimodal frequency 
distribution. Parameters are the same as in Fig. 1, except that now a = 0.4. (a) Analytical results 
for the evolution of |r(t)|; (b) numerical results; (c) comparison between both results; (d) evolution 
of the phase of the order parameter there is only a small time shift between analytical 

and numerical results. Notice that for an asymmetric frequency distribution with these parameter 
values, K > ‘Ija = 10/3, so that the synchronized phase is an asymmetric combination of clockwise 
and anticlockwise rotating traveling waves. 
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bution of Fig. 4 when wq = 200. The other parameter values are as in Fig. 4. 
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FIG. 6. Same as Fig. 4 for a lower value of the coupling constant, K = 4.5. Now 
Ka = 1.8 < 2 < 2.7 = K{1 — a). Only the component of the probability density with nega¬ 
tive frequency is synchronized. Then we obtain a traveling wave rotating clockwisely with constant 
|r(t)| and phase 'ip{t) = —ujot. 
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ip{t), for the unimodal Kuramoto-Sakaguchi model. Parameter values are: lvq = 15, /iq = 0.5, 
and K = 7.5. Notice the additional oscillation of the amplitude which is not predicted by our 
leading-order asymptotics. 
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FIG. 
that for 


9. Same as Fig. 7 for the following parameter values; loq = ho = 20, and K = 6. Notice 
such large values of the external field, a stationary state is reached for long times. 
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